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^ ■ Abstract 
O 

0^ . We study the convergence and the stability of fictitious dynamical methods 



for electrons. First, we show that a particular damped second-order dynamics 
has a much faster rate of convergence to the ground-state than first-order 



c . 

O . steepest descent algorithms while retaining their numerical cost per time step. 

Our damped dynamics has efficiency comparable to that of conjugate gradient 
methods in typical electronic minimization problems. Then, we analyse the 
factors that limit the size of the integration time step in approaches based 
on plane-wave expansions. The maximum allowed time step is dictated by 
the highest frequency components of the fictitious electronic dynamics. These 
can result either from the large wavevector components of the kinetic energy 
or from the small wavevector components of the Coulomb potential giving 
rise to the so called charge sloshing problem. We show how to eliminate 
large wavevector instabilities by adopting a preconditioning scheme that is 
implemented here for the first-time in the context of Car-Parrinello ab-initio 
molecular dynamics simulations of the ionic motion. We also show how to 
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solve the charge-sloshing problem when this is present. We substantiate our 

theoretical analysis with numerical tests on a number of different silicon and 

carbon systems having both insulating and metallic character. 
71.10+x, 71.20Ad 
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I. INTRODUCTION 

The introduction of a fictitious dynamics for the electron^!! with driving forces obtained 
from the total energy within density functional theory (DFT)I has provided a convenient 
approach to minimize the total energy of condensed matter systems and to perform ab-initio 
molecular dynamics simulations of the ionic motion. These techniques have been applied 
successfully to a variety of insulating, semiconducting, and metallic systems involving a large 
number of atoms in the context of structural optimization problems at zero temperature and 
of dynamical simulations of the atomic motion at finite temperature^. 

It is a subject of current interest to study the factors that limit the efficiency of fictitious 
dynamical methods for electrons in order to improve their numerical efficiency. This depends 
on the choice made for the dynamics and on the size of the time step that can be used to 
integrate numerically the equations of motion. 

When discussing how to choose a specific dynamics, it is convenient to consider total 
energy minimization separately from molecular dynamics. It has been shown by Car and 
Parrinello that to simulate the classical adiabatic motion of the atoms it is useful to adopt 
a second-order Newtonian dynamics also for the electronic degrees of freedom, since this ex- 
ploits optimally the concept of continuous simultaneous evolution of electronic and atomic 
degrees of freedomBi. Newtonian dynamics conserves energy. Different approaches should 
be used to minimize the electronic energy, as it is required to start a molecular dynamics 
simulation or to solve an optimization problem at zero temperature. The simplest approach 
to minimization is provided by steepest-descent dynamics, which can be viewed as a dy- 
namics of the first-order in the time-derivativei. Steepest-descent dynamics, which requires 
only knowledge of the gradients of the energy functional, is not very efficient particularly in 
metallic situations. Better schemes require some knowledge also of the second derivatives 
of the energy functional either explicitly or implicitly. Conjugate gradient methods have 
been developed in this contexti^ and have been shown to be superior to steepest-descent 
methods, particularly when the full energy functional was used in the line minimizations 
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and full account was taken of the orthonormality constraints on the wavefunctionsO. 

In this paper we show that a minor modification of a steepest-descent algorithm, namely 
replacing first-order dynamics with a specific damped second-order dynamics, improves sub- 
stantially the rate of convergence of the wavefunctions to the ground-state. The resulting 
scheme, which we call damped molecular dynamics, has efficiency comparable to that of the 
best conjugate gradient algorithms when used in typical electronic minimization problems, 
with the additional advantage of having basically the same numerical complexity of simple 
steepest-descent algorithms. 

We then investigate what determines the maximum allowed time step for numerical in- 
tegration when using steepest-descent (SD), damped (D) or Newtonian molecular dynamics 
(MD). In all cases the time step is limited by the need to integrate the high frequency com- 
ponents of the fictitious dynamics. These arise either from the large wavevector components 
of the electronic kinetic energy or from the small wavevector components of the Hartree 
energy due to the divergence of the Coulomb potential at small wavevector. In the latter 
case the related numerical instability is usually referred to as the 'charge sloshing' problem 
and it is expected to become serious when the size of the system becomes very large. 

The large wavevector instability can be eliminated by preconditioning the equations of 
motion since at large wavevectors the wavevefunctions are dominated by the kinetic energy 
and are to a large extent free-particle like. Indeed, it was already suggested earlier by 
several authors that this property could be used to speedup iterative schemes for electronic 
minimization. In particular, Ref. ( ^ proposed an analytical integration scheme for the large 
wavevector components of the wavefunctions within second-order dynamics. This scheme 
was subsequently extended in Ref. ( [1^) to first-order steepest-descent equations. Since this 
approach can be less stable than standard steepest-descent algorithm^ we will not discuss 
it any further. So far the most successful preconditioning scheme in the context of electronic 
minimization is the one proposed in Ref. ( |^) in connection with a conjugate- gradient method 
to minimize the total energy. 

In this paper, we propose a preconditioning scheme which is appropriate to all the dy- 
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namical methods referred to above, namely SD, D, and MD dynamics. It consists in properly 
scaling the fictitious masses associated to the large wavevector components of the electronic 
wavef unctions, in order to compress the high frequency spectrum of the electronic dynamics 
and to use a larger integration time step. Our preconditioning method is similar in spirit 
to the one described in Ref. ( 0) in the context of conjugate gradient minimization but it is 
formulated as a modification of the differential equations leading to SD, D, and MD dynam- 
ics. In particular, we apply it here for the first time to the Car-Parrinello MD equations, 
which provide an efficient approach for ab-initio molecular dynamics simulations of the ionic 
motion. In this context our preconditioning scheme allows to use a timestep which is two to 
three times larger than in previous applications of this method, resulting in a considerable 
saving of computational time. 

We now turn our attention to the 'charge sloshing' problem. This has been discussed 
previously in the context of self-consistent diagonalization of the Kohn-Sham Hamiltonian0. 
The onset of this kind of instability is expected to occur at significantly larger sizes in the 
context of fictitious dynamical methods since in these approaches the wavefunctions change 
little over a single timestep. Indeed recent MD simulations for metallic liquid silicon have 
shown no sign of a 'sloshing' instability up to cubic cells containing 216 silicon atoms'. 
However one expects that for sufficiently large cells the 'sloshing' instability should appear, 
although a quantitative theoretical analysis of it in the context of fictitious dynamical meth- 
ods for electrons has been so far missing. 'Sloshing' instabilities have been found numerically 
within some iterative schemes for electronic minimization in the case of systems having a 
very long linear dimension!!. In this paper we present a theoretical analysis of the 'charge 
sloshing' problem in the context of SD, D, and MD equations of motion. We find that the 
'sloshing' instability is absent for insulators, but it is present for metals. This is in accord 
with previous results of Ref. ( |Tl|) . A practical scheme to control the sloshing instability is 
discussed in the paper. 

To summarize, we improve the numerical efficiency of fictitious dynamical methods for 
electrons in several ways. Firstly, we replace steepest descent dynamics by a more efficient 



damped second order dynamics to minimize the total energy. Secondly, by preconditioning 
the fictitious electronic masses we increase the integration time step for total energy mini- 
mization and for simulation of the adiabatic ionic dynamics. Thirdly, we show that in the 
context of fictitious dynamical methods the so called 'charge sloshing' problem, which is 
expected to arise for large systems, is less serious than expected. We support our theoretical 
analysis with detailed numerical tests on several systems involving Si and C atoms. 

The paper is organized as follows. In Sec. II we discuss first-order SD dynamics and 
second-order conservative MD dynamics for the electronic degrees of freedom. In Sec. Ill 
we introduce a damped second-order dynamics which is substantially more efficient than 
SD and is competitive with the best conjugate gradient schemes for electron minimization. 
In Sec. IV we discuss large wavevector instabilities and the 'charge sloshing' problem. In 
Sec. V we discuss the preconditioning of large wavevector components. In Sec. VI we 
present some details of the numerical implementation. Finally, in Sec. VII we present the 
results of realistic numerical tests on silicon and carbon systems. Sec. VIII is devoted to 
our conclusions. 



II. FICTITIOUS DYNAMICS FOR THE ELECTRONS 

Dynamical methods for minimizing the electronic total energy and for simulating the 
adiabatic motion of the atoms are based on a fictitious dynamics of the electronic degrees 
of freedom. Within these approaches the forces acting on the electronic degrees of freedom 
are derived from the total electronic energy in the DFT-LDA form: 

E[{ij}] = Ek,n[W] + E,,t[p\ + ^h[p] + ^xc[p], (1) 

where Ekm, Eext, Eh, and E^c denote kinetic, external potential, Hartree and exchange- 
correlation energy, respectively!. The local density approximation is adopted for the latteii. 
The electronic charge density p is given by: 

p(r) = 2 < ^,|r >< r|^i >, (2) 



where the occupied orbitals j^j > are orthonormal. The factor of 2 accounts for the occupa- 
tion numbers, which here and in the following are supposed to be all equal to 2. Summation 
over repeated indices is understood. 

In order to ensure the orthonormalization of the electronic orbitals during a dynamical 
evolution it is convenient to add appropriate forces of constraints. These do not perform 
work on the electronic system, and can be conveniently calculated in terms of Lagrangian 
multipliers. The corresponding equations of motion for the first order dynamics are: 

Ml^<>=-i™ + A.,IV,>, (3) 

and those for the second order dynamics are: 

/.*>=-i™ + A«l*,>. (4) 

where we have assumed that the wavefunctions \ip > are real. Here and in the following the 
indices i and j run over the occupied states only. The symmetric matrix Ajj of the Lagrange 
multipliers enforces the orthonormality condition, i.e. < >— 5ij. The derivatives of 

£■[{■0}] with respect to the j^'i > define the Kohn-Sham Hamiltonian: 

The parameter // is a fictitious electronic mass. It is used to tune the speed of the electronic 
dynamics and does not describe any other physical property. When the ions are held fixed, 
this mass can be included in the definition of the time step and it is irrelevant. However when 
we allow the ions to move, the ratio between /i and the physical ionic masses is important 
since it defines the relative speed of the ionic and of the fictitious electronic motion. 

Now let us suppose that the wavefunctions are close to the minimum of the energy 

\A >= |# > >, (6) 

where {ipf^^ > are the wavefunctions at the minimum and {Sipi > are the corresponding 
deviations. We notice that \Sipi > have to fulfill the orthonormality condition to linear 
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order, i.e. < ipl l^ipj >= 0. To linear order in \Sipi >, the Lagrange multipliers Aij are the 
same for first and second order dynamics, and are given by 

Aij = 2 < ^pilHKsli^j > ■ (7) 

Thus, by retaining only the terms up to linear order in {Sipi > in the equations of motion 
(I) and (^) we obtain: 

\6ij, >= -k,,\5ip, > (8) 

and 

\5^i >= -Ki,\Sijj >, (9) 

respectively. Here Kij is a linear operator, which acts on the single-particle Hilbert space 
and which has the same form for both first- and second-order dynamics. In the following we 
will use the notation K to indicate a matrix of operators having for elements the Kij. Notice 
that K is a positive definite linear operator since Eqs. (H) and (|^) result from a quadratic 
expansion of in \6tpi > about the minimum i?[{'?/'o}]) i-e.: 

EiW] - E[{^o}] =f^< StP,\k,^\6^^ > +0{6i;'). (10) 

The equations of motion (H) and can be formally integrated yielding: 

\5ij^{t) >= (expi-Kt)) . . |5^,(0) > (11) 

and 

|5^.(t) >= {cos^t)ij\64jj{0) > +{K~^/'^sin\fkt)ij\6ijj{0) > , (12) 

respectively. In the case of first order evolution the wavefunctions decay exponentially to- 
wards the minimum -E[V'o]; while in the second order evolution they perform small oscillations 
around it. These motions take place with characteristic decay rates and frequencies which 
are equal to the eigenvalues Ka of the operator K and to the square root y/K^ of these 
eigenvalues for first and second order dynamics, respectively. 
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In numerical implementations the electronic states are expanded on a finite basis set, so 
that only a finite number of eigenfrequencies and eigenmodes occur. Let Kmin and K^ax 
be the minimum and maximum eigenvalues of K. The maximum allowed time-step for 
numerical integration is proportional to the smallest period of the system, i.e. to l/K^ax or 
to 1 / a/ Kmax for first and second order dynamics, respectively. 

In the case of first order dynamics, the minimum eigenvalue dominates the long-time 
behavior of the decay to the ground-state, so that a rough estimate of the convergence time 
is given by 1/Kmin- Recalling that the size of the time step is proportional to 1/Kmax, one 
finds that the number uqi of integration steps needed to converge satisfies the condition: 



In the case of second order dynamics we usually start a simulation from an electronic 
configuration close to the minimum of the electronic energy. Then if the ionic and the 
electronic frequencies are well decoupled0i, the electrons remain adiabatically close to the 
instantaneous energy minimum during the ionic evolution. Let Uion be a typical ionic fre- 
quency. The adiabatic condition requires it to be much smaller than the minimum electronic 
frequency i.e. 



A meaningful measure of the simulation's workload is given by the number of time steps 
necessary to integrate a full ionic oscillation. Thus, recalling that the time step is inversely 
proportional to y/Kmax, we find that the number no2 of steps necessary to integrate a typical 
ionic oscillation satisfies the following condition: 



noi oc 



(13) 




(14) 




(15) 



III. DAMPED SECOND ORDER DYNAMICS FOR MINIMIZATION 



In the previous section we showed that the number of iterations necessary to minimize 
the electronic energy within steepest descent dynamics is proportional to K^ax/ K^in. In 
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this section we present an improved minimization dynamics in which the typical number of 



iterations is instead proportional to \J K^ax/ Kmin, i-e. a number significantly smaller than 
Kmax/ Kmin- We attain this goal by inserting in Eq. (|D a damping term as follows: 

>= -^^^^ - 27/^1^^ > +^ij\^J > (16) 

This equation defines a damped second order dynamics. As in the previous section we 
study the resulting motion close to the energy minimum. We find that the deviations of the 
wavefunctions from the minimum are subject to damped oscillations given by: 

ISipiit) >= exp(iW+t)ij#]+^ > +exp(iW_t)ij#j"^ > (17) 

Here iV'j^'* > s-nd iV'j "* > are determined by the initial conditions, and 



W± =27l± VK-721 (18) 

The real part of W± gives the frequencies of the oscillatory motion, while its imaginary part 
gives the decay rate to the minimum. In order to maximize the rate of convergence, we must 
use the maximum value of 7 for which the argument of the square root remains positive. 
This optimal value of 7 is given by 



7opt = yKmin, (19) 

since this value corresponds to critical damping of the smallest eigenvalue of K. In this case 
the imaginary part of all the eigenvalues of Wj. is equal to 'jopt and the time of convergence 
to the minimum is of the order of Xj \/Km%n- The integration time step is related to the 
maximum norm of the eigenvalues of Wj., which is equal to s/K^ax- Thus, the number of 
integration steps necessary for minimization is given by: 



nD02 oc J oc ,Jn^. (20) 

V ^iriin 

From this formula we see that a relevant gain of efficiency is obtained when using damped 
dynamics instead than steepest-descent dynamics to minimize the electronic energy. The 
gain is particularly important when a large number of iterations is needed to converge to 
the ground-state, which is typically the case of metallic systems. 
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IV. SECOND ORDER EXPANSION OF THE LDA ENERGY FUNCTIONAL 



In this section we compute explicitly the eigenvalues of the operator K. For this purpose 
we consider the expansion of the energy functional around its minimum '^^'^^ up to second 
order in 6"$. This is given by: 

-^[{^W}] = 2 < 6^,\HKsm^ > -2 < 6^,\6^, >< ^f^lHKsl'^f > (21) 



+ J dr J dr'5p{r) 



+ 



5p(r') + 0({5^3}). 



6p{r)5p{r') 6p{r)5p{r') 
The second term on the r.h.s. of this equation comes from the Lagrange multipliers (see Eq. 
(0) ), and 

6p{r) = 4 < V >< r|5^,: > (22) 

gives the variation of the electronic density to first order in {5^E'}. We recall that {\E'} and 
are supposed to be real. 

A. Non self-consistent case 

If we neglect the last two terms in Eq. (pT]), i.e. the terms corresponding to variations 
of the Hartree and exchange-correlation potentials, we recover the expansion of the total 
energy appropriate to a non-selfconsistent Hamiltonian. Then we can expand the ^^^^ and 

in terms of the the real eigenvectors Ixl' > of Hks^ which have eigenvalues Since the 
total energy is invariant under unitary transformations in the subspace of occupied states, 
we can suppose without loss of generality: 

l^f >= |X° >, 
|5^.>=E4lx2>, (23) 

k 

where are real coefficients. Here and in the following the indices i and k refer to occupied 
and unoccupied states, respectively. Hence, as shown in Ref. ( |13D, we obtain for 5Ensci i-e. 
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the second-order variation of the energy in which the selfconsistency of the potential is not 
taken in account: 

6Ensc = 2J2{c'k?{ek-e^). (24) 

ik 

where £j and Ek are respectively the occupied and the unoccupied eigenvalues of H^s- By 
comparing Eq. (|2^) with Eq. ([l^) we see that the eigenvalues of K are given by 

= 2^^^ (25) 

and the lowest eigenvalue of K is given by Kmm = '^Egap/ fi, in terms of the energy gap Egap 
separating the lowest unoccupied from the highest occupied electronic level. 

In the case of an insulator, the energy gap has a finite positive value which, above a 
certain size, is independent of the simulation cell. In the case of a metal instead, the energy 
gap is still finite and positive for a finite sized system but it is no longer independent of 
the simulation cell. In fact the energy gap and K^m tend to zero for a cell size going to 
infinity. However, many properties of interest do not require an infinite energy resolution 
for the states around the Fermi energy. Typically a small but finite energy resolution Eerr 
is sufficient. E^rr does not depend on the size of the system, and K^rr = '^E^rr/fJ' replaces 
Kmin in Eq.s (p!3l),(p0|) and (|^) to estimate the convergence rates uqi, nD02 and the optimal 
damping parameter •jopt- Since E^rr is much smaller than a typical energy gap of an insulator, 
the number of iterations needed for ground state convergence is much larger in metals than in 
insulators. For the same reason, a perfectly adiabatic separation between ionic and electronic 
motions is not possible for metals. However, as shown in Ref. ( [1^) a satisfactory solution 
to this problem, in the context of Car-Parrinello simulations, can be obtained by using two 
Nose' thermostats to control separately the respective temperatures of the ions and of the 
electrons. 

When expanding the wavefunctions in terms of plane-waves, we can define an effective 
cut-off energy E^ut given (in a.u.) by q^axl'^^ where q^ax is the largest wavevector in the 
basis set. The band of empty states is usually much larger than the band of occupied states. 
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Thus when in Eq. (^) the index k refers to the highest unoccupied states, the eigenvalues 
/^(j of K have a neghgible dependence on the occupied state index i. Furthermore, since 
the highest unoccupied states are free-particle like, the energy difference Ek — Si is dominated 
by the kinetic energy of the state k, i.e. 

Sk-e^^ gV2, (26) 

where q is the wavevector associated to the state k. The maximum eigenvalue of K is 
therefore approximately given by Kmax — '^{imax/'^) / 1^ = 2i?cut/At- It is this eigenvalue that 
limits the maximum allowed time step for numerical integration in a non self-consistent case: 
the numerical integration becomes unstable and the time step has to go to zero when E^ut 
goes to infinity. 



B. Self-consistent case: charge sloshing 



We now consider the terms of Eq. (|2T| ) that we neglected in the previous subsection in 
order to see if they affect the maximum eigenvalue of K. In this case K is not diagonal in 
the representation of the c\ and, for an arbitrary system, it is not possible to diagonalize it 
analytically. Thus we need some simplifying assumptions. Let us consider a crystal of given 
periodicity and use a supercell containing an arbitrary number of replicas of the crystal 
unit cell. In this case the {^E'o} are linear combinations of Bloch-functions with the crystal 
periodicity whereas the fluctuations {5'^} may have all the wavelengths compatible with the 
supercell. In other words we are restricting the periodicity of the unperturbed state but not 
the periodicity of the fluctuations. Based on the above simplifying assumption we find that 
charge sloshing affects differently metallic and non-metallic systems. A numerical example 
presented in Sec. VII suggests that this result should hold also for non-periodic systems. 

In order to find out whether the maximum eigenvalue of K diverges when the supercell 
size tends to infinity, we restrict our analysis to the Hartree term since the LDA exchange- 
correlation energy is well behaved and typically has a negligible effect compared to the 
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kinetic energy on the maximum eigenvalue of K. The second order variation of the Hartree 
energy is given by: 

6Eh = jdvj rfr'5p(r)--L-5p(r') = E E ^^7^^I'^P(P + (27) 

I I G- P I I 

Here p is a vector belonging to the first Brillouin Zone of the crystal, G is a vector of 
the reciprocal lattice of the crystal and the sums extend over all the non-zero wavevectors 
p + G = q — q' where q and q' are two generic plane- waves of the basis set used to represent 
the electron wavefunctions in the supercell of volume VL. 5p(p + G) is the Fourier Transform 
(FT) of 5p{r). 

When a linear dimension L of the supercell becomes very large, (p + G)^^^ = {p)min, i-e. 
the smallest nonzero q — q' vector, tends to zero like 1/L. If, correspondingly, the maximum 
eigenvalue of K diverges, we have the so-called 'charge sloshing' scenario. To study the 
effect of Pmm on the maximum eigenvalue of K we consider only the terms with G = in 
Eq. (pTl). Then since Sp{r) is real, |5p(p)| = |5p(— p)| and 6Eh{G = 0) can be written as: 

6Eh{G = 0)= E 7?ilMp)r, (28) 

where 

6pip) = E4(< x^\cosipr)\xl > 4 + ^ < x^\stnipv)\xl > 4) (29) 

ik 

Notice that Eq. (^Sj) is a quadratic form in terms of the 4- The real coefficients 4 can 
be considered as the components of a real vector |c >. Similarly we can introduce two 
real vectors |A(p) > and |B(p) > whose components, labelled by the composite index (ik), 
are given by < Xi I cos (pr)|x° > and by < Xi k^^(pi')|Xfc >; respectively. In this notation 
Sp{p) = 4(< c|A(p) > + i< c|B(p) >) and Eq. (|2|) becomes: 

6Eh{G = 0)= E |^(<c|A(p)>< A(p)|c> + <c|B(p)><B(p)|c>). (30) 

Using the fact that the are eigenstates of a periodic crystal, it is easy to showlil that the 
vectors |A(p) > and |B(p) > constitute an orthogonal set: 
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<A(p)|A(p')> = Wf^(p) 

<B(p)|B(p')> = Wf5(p) (31) 
<B(p)|A(p')> = 
where PxiV'x > 'static structure factor' S'(p) is defined by: 

^(P) = < x'Ae-'nxl X Xlle'-nx'^ > • (32) 

ik 

Hence |A(p) > and |B(p) > are the vectors that diagonahze the quadratic form in Eq. (pO|). 
The corresponding eigenvalues of K are given by: 



1287r ^, , 1287r ^, , GAir S{p) 



Therefore, when a hnear dimension L of the supercell tends to infinity and, correspondingly, 
p^i„ goes to zero, i^A(p^„) and i^B(p„,j do not diverge if S{p) is of order Oip"^). 

We now consider a jellium model as a representative metallic system. In this case the 
x'i are plane waves and one finds0: 5'(p) = p[l — {p/pF)'^/l'2\pp/8rr'^, where pp is the Fermi 
momentum, and p < 2pp. As a consequence for L going to infinity, _ft'A(p,„-,j and -ft'B(p„.„) 
diverge as L and the time step for numerical integration has to be reduced accordingly: this 
is a charge sloshing situation. 

When the system is a periodic insulator one finds instead that S'(p) goes to zero as p"^ 
(see Appendix). As a consequence for L going to infinity, -^A(p„-,j and Kb{p^-j tend to a 
constant and the time step for numerical integration is independent of L: charge sloshing is 
absent here. 

We stress that the above conclusions apply only if we consider small fluctuations around 
the ground-state: this is the typical case of ab-initio molecular dynamics simulations of 
the ionic motion. However, in the initial steps of an electronic minimization procedure, 
the wavefunctions may be far from the ground-state. In this case it is possible to observe 
sloshing effects also in periodic systems having an insulating ground-state. 

Since charge-sloshing is a consequence of the singularity of the Coulomb potential at 
small p, a simple way of eliminating charge-sloshing instabilities consists in replacing the 
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Coulomb potential Air/p'^ with a Yukawa potential 47r/(j9^ + a^), where 27i'/a is a typical 
decay length of the order of the system size that corresponds to the onset of the sloshing 
instability. In the case of an insulator we can use this technique to stabilize the numerical 
integration during the initial steps of an electronic minimization run. Then when we are 
sufficiently close to the ground-state we can set a = and converge to the exact ground- 
state. We will show with a numerical example in a subsequent section that this technique 
allows us to converge to the exact ground-state of a disordered insulating system with a 
number of steps independent of the system size. In the case of a metal it is not possible to 
set a equal to zero not even in the proximity of the ground-state. However we notice that 
Lmin is usually much larger than the typical screening length of a metal. The results of a 
numerical simulation for a large but finite metallic system should not change appreciably 
if the Coulomb potential is replaced by a Yukawa potential that is equal to the Coulomb 
potential for distances smaller than Lmin- 



V. PRECONDITIONING THE EQUATIONS OF MOTION 

The numerical efficiency of all the fictitious dynamical methods previously introduced 
can be improved by preconditioning the dynamics in order to reduce the ratio Kmax/Kmin- 
This can be achieved by replacing the constant fictitious mass parameter /i in Eqs. (^^, |16|). 
with an arbitrary positive definite operator fi. The resulting increased arbitrariness in the 
choice of fi can be exploited to compress the highest frequency components of the spectrum 
of the fictitious electron dynamics. Recalling that these are due basically to the high energy 
unoccupied states which are free-particle like (see Eq. (P^D), we choose an operator fi which 
is diagonal in g-space with eigenvalues given by: 

= /io if < Ep ^^^^ 

Below a certain cutoff energy Ep, it is worth considering a constant mass /io, because the low 
energy eigenstates have a relevant potential energy contribution and are not free-particle like. 
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The preconditioning cut-off Ep therefore represents the threshold above which the states are 
dominated by the kinetic energy. 

It is easy to show that the solutions of the preconditioned equations of motions for small 
displacements, are still given by Eqs. (|Tl],|T2|,[l^) if the operator K is replaced by the operator 
K characteristic of the preconditioned dynamics. All the relations ([TB|), (|TB|), (^Dj) found for 
first and second order dynamics with and without damping hold therefore also in the case 
of the preconditioned dynamics but, in the latter case, Kmax and Kmm have to be replaced 
by the maximum and minimum eigenvalues of K, i.e. by Kmax and Kmin- 

The preconditioning cut-off Ep that minimizes the ratio Kmax/ is called the optimal 
preconditioning cutoff. It depends strongly on the atomic species, i.e. on the pseudopoten- 
tials and on the plane-wave cutoff that are used in the calculation. It depends only negligibly 
on the physical environment. Thus, for a given atomic species, it is possible to find the op- 
timal preconditioning cutoff by performing calculations on a simple reference system. We 
present a typical example in Sec. VII. 



VI. NUMERICAL IMPLEMENTATION 



In our numerical implementation we adopt the standard procedures described in Refs. 
to integrate the equations of motion for first and second order dynamics. In the 



case of damped second order dynamics we follow the procedure introduced in Ref. ( to 
integrate Car-Parrinello dynamics in presence of a Nose- Hoover thermostat. We obtain for 
first, second order and damped dynamics, respectively: 

|^i(t + A) >= \iPi{t) > -2js-^HKs\^i > A + Xi,-/i-Vj W > (35) 

\tljiit + A) >= - A) > +2|^i(t) > -2fi-^HKs\^i > A^ + X,.jfi-^\i,.j{t) > (36) 

\^i{t + A) >= > + 

+ {\Mt) > -mt-^) > -fi-'HKs\A > ^) ^-^^+X,,ix-'\^jit) > (37) 

17 



where A is the integration time step and Xij is a symmetric matrix equal to AyA for 
first order dynamics, equal to Aj^-A^ for conservative second order dynamics, and equal to 
AjjA^/(l + 7A) for damped second order dynamics. The matrix X is found by imposing 
the orthonormality of the wavefunctions at time t + A: 

<^i{t + A)\^j{t + A) >= 5ij (38) 

We notice that the inversion of the mass operator is straightforward in q space where it is 
diagonal. For the calculation of X^j we define the wavefunctions \ijji{t + A) > as the r.h.s. 

without the orthonormalization terms Xijfi~^\iljj{t) >. Then Eq. ( |38D 
becomes: 

XMX^ + BX^ + XB^ = 1 - A (39) 

where the matrices M, B, A are given respectively by: 

Mi, =< i;i{t)\il-^\ijjit) > (40) 



oftheEq. (B5,3|,37 



Bi,=<^i{t + A)\fi-'\tljj{t)> (41) 

Aij =< i,i{t + A)\i,j{t + A) > . (42) 

The scalar products are easily evaluated in g-space where the mass operator (i is diagonal. 

Eq. (|39|) is formally identical to the matrix equation that expresses the orthonormality 
condition for Car-Parrinello dynamics when using Vanderbilt's ultrasoft pseudopotentials0. 
It can be solved as described in Ref. ( [1^). The matrix B can be conveniently split into a 
symmetric part Bg and an antisymmetric part Ba. The antisymmetric part Ba is first order 
in A, while X and 1 — A are first (second) order in A for first (second) order dynamics. 
Using these properties, we can solve Eq. ( pOf ) iteratively in terms of increasing powers of A 
(A^): 

B.X^'^+i) + X("+i)Bs = 1 - a - X(")MX(") - BaX(") + X(")Ba (43) 



Here X*^*^) is the solution of the equation: 

BsX(°) + X(°)Bs = 1 - a (44) 
and the l.h.s of the Eq. (^ , ^4|) is inverted after transforming to the basis where Bg is 



diagonal 

An alternative approach based on unconstrained total energy functional which avoids 
explicit orthonormalization has been recently proposed in Refs. ( 0,^). The electronic 
mass preconditioning scheme discussed in the present paper can be easily applied to the 
unconstrained energy functional method without any overload. 

In some applications of fictitious dynamical methods for electrons the orthonormaliza- 
tion of the electronic wavefunctions can be achieved via a Gram-Schmidt procedure. We 
stress that this approach is not justified in connection with the mass preconditioning scheme 
described above. Indeed if the Gram-Schmidt orthogonalization procedure is used within 
preconditioned steepest-descent dynamics, one is not guaranteed that the energy will de- 
crease at any integration step for a sufficiently small time step. The origin of this instability 
is related to the non-holonomic character of the constraints imposed via a Gram-Schmidt 
procedure. We found that this instability is rather severe in practical numerical applications, 
where it spoils all the efficiency gains of the mass preconditioning scheme. 

VII. NUMERICAL RESULTS 

We tested the different dynamical schemes described above on various physical systems 
within a DFT-LDA formulation. In particular we considered Si and C systems. We used 
pseudopotentials of the Bachelet-Hamann-Schliiter typeS, with s and p non-locality in the 
Kleinmann-Bylander form0. The cut-off for the plane-wave expansion of the electronic 
orbitals was 12 Ry for silicon, and 35 Ry for carbon. We carried out all the calculations 
at the r point of the Brillouin zone only. Moreover, in order to compare the dynamical 
schemes with conjugate gradient minimization, we used a tight-binding energy functional 
for carbon@. 
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A. preconditioning 



We start by presenting the results obtained with preconditioning. In order to deter- 
mine the optimal preconditioning cut-off Ep we had to minimize the ratio Kmax/ Kmm- We 
measured Kmax and Kmin within first order and second order dynamics, by giving a small 
displacement to the system from its energy minimum. In the case of first order dynamics, 
the numerical integration of Eq. (P^D, becomes unstable and results in an exponential in- 
crease of the energy, when A > 2/ Kmax- Therefore the maximum allowed integration time 
step provides an accurate way of estimating K^ax- Kmin, i-e. the lowest eigenvalue of K, 
gives instead the slowest rate of decay of the energy. This rate is conveniently sampled at 
large times, i.e. when only the slowest exponential is left in the decay. 

We report in Fig. (1) the ratio Kmax/ Kmin as a function of Ep for a Sis molecule. We 
notice that for the highest Ep values the ratio decreases linearly with decreasing Ep. The 
behavior of the ratio Kmax/ Kmin in this range is explained by the following considerations. 
First, the minimum frequency is unchanged, since it is related to the lowest excited state 
which has small components at high q. Second, all the excited modes at energies higher than 
Ep are compressed to the same maximum frequency 2Ep/ fiQ, as long as they are kinetically 
dominated. Instead, in the range of low preconditioning cut-offs Ep, the minimum frequency 
Kmin decreases and the highest excited modes become less efficiently compressed. Thus a 
minimum value of the ratio Kmax /Kmin is found, as we can see in Fig. (1). This minimum 
occurs at Ep=l Ry. The corresponding reduction of the ratio Kmax/ Kmin is of a factor of 5 
compared to the non-preconditioned case. We obtained very similar results for a sample of 
crystalline silicon in the diamond structure, where the optimal Ep was also close to 1 Ry. 

In the case of second order dynamics \J Kmax and \J Kmin can be found as the maximum 
and minimum frequencies of the power spectrum of the fictitious electronic dynamics. This 
is easily evaluated by computing the velocity autocorrelation function corresponding to 
the wavefunction dynamics. The power spectrum of the velocity autocorrelation function 
corresponding to the electronic fictitious dynamics is given in Fig. (2) for the case of the 
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Sis molecule. In particular we show results obtained with optimal preconditioning {Ep = 1 
Ry) and unpreconditioned dynamics. In calculating the spectra we chose the value of fiQ in 



such a way that the lowest frequency J Kmin of the preconditioned dynamics coincided with 



\/Kmin, i-e. the lowest frequency of the dynamics without preconditioning. This is achieved 
by setting /io=260 a.u. when the mass associated to the dynamics without preconditioning 
is /i=300 a.u.. The significant compression of the high frequency modes resulting from 
preconditioning is clearly evident in Fig. (2). 

For a sample of crystalline carbon in the diamond structure we found an optimal value 
of Ep equal to 2.7 Ry. This reduced by a factor of 9 the ratio Kmax/Kmin compared to 



the unpreconditioned case. In this case in order to make the lowest frequencies J Kmin and 



\JKmin to coincide, the mass /i of the unpreconditioned dynamics had to be rescaled by a 
factor of 0.93 in order to obtain the mass /iq of the preconditioned dynamics. 

We notice that in a different context the authors of Ref. ( |^) proposed to use a precondi- 
tioning cutoff Ep equal to the expectation value of the kinetic energy divided by the number 
of electrons. In the cases discussed above this corresponds to a value of Ep = 0.8 Ry and 
Ep =2 Ry for Si and C, respectively. These values are close to the optimal values of Ep. 



In order to test the effect of preconditioning on ab-initio molecular dynamics simulations 
of the ionic motion, we considered the coupled set of equations given by Eq. (H) for the 
electronic degrees of freedom and by 



for the ionic coordinates Rj. Here Mj are the physical ionic masses and the mass fi in Eq. 
(^) has to be replaced by the mass operator fi in the preconditioned case. Eq. (|^) and ( ^Sf ) 
reproduce the adiabatic dynamics of the ions when the appropriate decoupling condition, 





B. Ionic Molecular Dynamics 



(45) 




II, is satisfiedlirH. We considered the vibrational motion of 




a 
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Sis molecule during a time span of about 0.3 ps. In the unpreconditioned case we used a 
time step A=7 a.u. to integrate the equations of motion. This is close to the maximum 
allowed time step for a fictitious electronic mass /i=300 a.u.. Preconditioning allowed to 
increase this time step to A =15 a.u. for a mass yUo=260 a.u. and a preconditioning cut- 
off Ep=l Ry. In spite of the significantly larger time step the preconditioned dynamics 
proceeded adiabatically in the same way as the one without precontioning. In particular, 
any systematic energy transfer from the ionic system to the electronic one was absent. We 
plot in Fig. (3) the temporal evolution of the ionic kinetic energy and of the longest side of 
the Sis molecule as a function of time in both the preconditioned and the unpreconditioned 
cases. Differences between the two dynamics are not noticeable. 



C. Damped dynamics in Insulators 

In order to assess the efficiency of the various minimization dynamics discussed in this 
work, we considered a 64 atom amorphous Si sample generated by ab-initio molecular 
dynamics^. We notice that this system has a finite gap, and therefore a non-zero Kmin- 
In all our total energy minimizations we used the same set of starting trial wavef unctions. 
These were obtained by minimizing the total energy with a very small energy cut-off Ecut 
of 2 Ry. We then minimized the total energy with a cut-off of 12 Ry using four types of 
dynamics, namely steepest descent and second order damped dynamics both without and 
with optimal preconditioning. We report the results in Fig. (4). In particular, we found 
that, when using the optimal value 'jopt of Eq. (p!9D , the rate of convergence of second order 
damped dynamics is faster than that of steepest descent dynamics by the amount expected 
from the theoretical analysis in Sec. III. Preconditioning accelerated further the rate of 
convergence, so that finally the rate of convergence of preconditioned second order damped 
dynamics was 14 times faster than the one of unpreconditioned steepest descent. 

We determined the value •jopt by a rough estimate of Kmin based on steepest descent 
dynamics. In particular, a 3 point fit of the exponential decay of the total energy gives: 
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where Ei E2 E3 are the energies at three successive steps of steepest descent. We waited 
until only the slowest exponential was left in the decay. If faster exponential are still present, 
Eq. ( ^61) overestimates 7optA. In a practical calculation, we therefore suggest to start the 
minimization with a few steps of steepest descent and to use Eq. ( ^6]) to obtain an upper 
bound for the optimal 7 A. Then, we suggest to proceed with the damped second order 
minimization, readjusting 7A in order to achieve the optimal limit of critical damping. 
As we can see from Fig. (4), it is indeed convenient to use steepest descent in the first 
steps of minimization when the highest frequency components dominate the deviation of 
the energy from the minimum. Subsequently, when only the slowest frequencies are left, 
damped dynamics becomes much more convenient, especially in those cases of utterly slow 
convergence rate. 



D. Damped dynamics in Metallic Systems 

In order to test the efficiency of our damped dynamical scheme for minimization in the 
case of metallic systems we applied it to liquid silicon which is a metal. We use a 64 atom 
sample generated by ab-initio molecular dynamicsEl. As explained in Sec. IV, the damping 
constant 7 can be fixed on the basis of the required energy resolution E^rr, for which we 
chose here a value of about 20meV. We minimized the total energy with a cutoff of 12 
Ry using four types of dynamics, similarly to what we did in the insulating case. Again 
the starting trial wavefunctions were obtained by a minimization using a small cutoff of 2 
Ry. We report the results in Fig. (5). Notice that in the present metallic case, steepest 
descent dynamics is particularly inefficient, while damped dynamics is very effective, since it 
improves by many orders of magnitude the convergence rate of steepest descent. A further 
gain results from preconditioning. 
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E. Comparison with Conjugate Gradient Minimization 



In this subsection we compare our damped dynamical method with a conjugate gradient 
minimization scheme. The standard conjugate gradient procedure, described e.g. in Ref. ( 
cannot be directly applied to a constrained functional, unless some additional simplifying 
assumptions are invoked which can reduce the minimization efficiencyi'i. To fully exploit 
the power of the conjugate gradient procedure the authors of Ref. ( ^ proposed to use an 
unconstrained energy functional. We adopt the same procedure of Ref. ( ^ but we use 
a different form for the unconstrained energy functional. We use the form suggested in 
in the context of electronic structure calculations with linear size scaling 
without imposing any localization constraints on the electronic orbitalJ^i. For reasons of 
numerical simplicity we adopt here a total energy functional based on a non self-consistent 
tight-binding Hamiltonian. This choice simplifies considerably the line minimization in the 
conjugate gradient scheme, which can be performed exactlyil. 

We used a tight-binding Hamiltonian for carbonS, and we considered an ionic liquid 
configuration of 64 atoms at a temperature of 5000 K. For this configuration the system is 
metallic. Our results are reported in Fig. (6) where we plot the logarithmic error in the 
total energy per atom versus the number of iterations for various minimization schemes, 
namely damped dynamics, conjugate gradient and steepest descent minimization. In the 
case of conjugate gradient minimization the number of iterations was multiplied by a factor 
of 2 in order to take into account the increase in computational cost arising from line 
minimization. From Fig. (6) it is evident that the numerical efficiency of both conjugate 
gradient and of damped molecular dynamics is considerably superior to that of steepest 
descent minimization. In the present example the numerical efficiency of conjugate gradient 
and that of damped molecular dynamics minimization are practically the same. 

We expect that the results that we have found here should remain valid also of in the 
case of a self-consistent LDA Hamiltonian. 



Refs. 
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F. Charge sloshing on very long cells 



In order to study charge sloshing effects we used a tetragonal supercell having a long side. 
In particular we considered crystalline silicon in the diamond structure and we constructed 
two supercells by repeating four or eight elementary cubic cells along the crystallographic 
(100) direction. The resulting supercells contain 32 and 64 atoms respectively. Then we 
broke the translational invariance of the diamond lattice by giving the silicon atoms a random 
displacement of about 5 percent of the bond length. This does not modify the insulating 
character of the system. 

In the present example we have considered only preconditioned steepest descent mini- 
mization. As in the previous subsections we prepared the initial trial state by minimizing 
the total energy with a small cutoff of 2 Ry starting from a set of random wavef unctions. 
Severe charge sloshing instabilities immediately showed up during this initial minimization 
in which the starting random wavefunctions were very far from the converged insulating 
ground-state. In particular, already for the 32 atom cell the time step for numerical integra- 
tion had to be reduced by an order of magnitude compared to the time step that we could 
use in an equivalent situation with a smaller cell. Such instability was completely eliminated 
by replacing the Coulomb by a Yukawa potential as described in Sec. IV. We adopted here a 
parameter I'nja = 20.5 a.u. for the Yukawa potential. Once obtained the initial trial state, 
we performed a total energy minimization on the 32 and on the 64 atom cell with a cutoff 
of 12 Ry. The results are shown in Fig. (7) which reports the deviation from the converged 
ground-state energy as a function of the number of numerical time steps. During the initial 
30 steps we used the Yukawa potential. This allowed us to use for both 32 and 64 atom cells 
an integration time step equal to the one usually adopted for the same system when using 
cells sufficiently small that no charge sloshing effects are present. Then we switched from 
Yukawa to Coulomb potential. After 30 minimization steps with the Yukawa potential the 
system was insulating and already very close to its exact ground-state. In this shown 
analytically in Sec. IV charge sloshing instabilities are not expected to occur in a periodic 
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system. Indeed not even in our disordered sample they did occur. Therefore we could use 
in the final 30 minimization steps the same time step used with the Yukawa potential. The 
overall convergence rate, as it can be clearly seen in the figure, is independent from the cell 
size. 

VIII. CONCLUSIONS 

We have presented a detailed analysis of the stability and of the convergence rate of 
fictitious dynamical methods for electrons. We have succeeded in improving considerably 
the efficiency of currently used algorithms for total energy minimization and for ab-initio 
molecular dynamics. 

In the case of ab-initio molecular dynamics simulations of the ionic motion we have 
introduced a novel preconditioning scheme which gives rise to an overall saving of CPU 
time of the order of 2-3 in typical applications. In the case of total energy minimization 
we have introduced an optimal damped preconditioned dynamics which has a convergence 
rate substantially faster than steepest descent algorithms and comparable to that of the best 
conjugate gradient schemes for electronic structure calculations. This is especially important 
in metallic situations. 

Although in this paper we confine our analysis to electronic minimization, we stress that 
the damped dynamics algorithm can also be applied to ionic minimization. In this case the 
optimal ionic damping parameter is related to the phonon frequencies of the system under 
study. 

In addition, we have presented a detailed analysis of the charge sloshing instability and 
we have indicated a practical way to control it. We have shown with a numerical example 
that, in the case of insulators, our approach allows us to converge to the ground-state with 
a number of iterations that is independent of the system size. 
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IX. APPENDIX 

In this appendix we show that for a periodic insulator the function S(p) given in Eq. 
(0) goes to zero like for p going to zero. 

S{v) = ^ < X°|e-*P^|x°fc X Xl\e+'nx'r > ^^^^ 

where the indices i and j refer to occupied states and the index k refers to empty states. 
In Eq. ( ^7|) we used the completeness relation: J2k \Xk >< Xk\ = 1 ~ I]j IXj >< Xjl- Let 
us suppose for simphcity that we have a single occupied band. Since the expression in Eq. 
(^) is invariant under unitary transformations on the occupied subspace, we can write it in 
terms of Wannier functions, i.e.: 

^(P) = 7T^(1 - E I < Wo\e-'nWn > n, (48) 

where is the Wannier function centered on site R, and fimm is the volume of the ele- 
mentary cell. The Wannier functions are exponentially localized in the case of an insulator: 
this allows us to expand in a Taylor series for p going to zero the exponentials in Eq. (^8|). 
In particular if we consider the term with R = in Eq. (|48|), and expand the exponentials 
in pr around p < r >= p < W^o|i"|Wo >, we get: 

1 - I < Wo\e-'P''\Wo > P = 

= 1 - I < Wo\l - tp{r- < r >) - [p(r- < r >)]V2|Wo > ^ + o{p'^) (49) 

= + < W^o|[p(r- < r >)]^\Wo > +o{p^) 
This term tends to zero as p"^. In a similar way one can show that the terms with R different 
from zero in Eq. (^Sf ) also go to zero as p'^. 
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FIGURES 



Kmax/ K^in ^ ^ function of the preconditioning cut-off Ep for the Sis molecule. A 
periodically repeated cubic cell of 20 a.u. is used in all the calculations for the Sis molecule. 

Spectra of the electronic frequencies for the Sis molecule. The solid hne refers to 2nd 
order dynamics without preconditioning. The dashed line refers to 2nd order dynamics with 
preconditioning {Ep—1 Ry). 

Ionic dynamics of the Sis molecule without (solid line) and with (dots) preconditioning. 
In (a) we report the oscillations of the long side of the Sis molecule, and in (b) the oscillations 
of the ionic kinetic energy as a function of time. 

Total energy minimization for a 64 atom amorphous silicon sample, using non- 
preconditioned steepest descent (SD NP), preconditioned steepest descent (SD P), non- 
preconditioned damped dynamics (D NP) and preconditioned damped dynamics (D P). We 
plot the logarithm of the difference between the energy per atom [E) and the ground state 
energy per atom {Eq) in Hartree units vs. the number of integration steps. 

Total energy minimization for a 64 atom liquid silicon sample, using non-preconditioned 
steepest descent (SD NP), preconditioned steepest descent (SD P), non-preconditioned 
damped dynamics (D NP) and preconditioned damped dynamics (D P). We plot the loga- 
rithm of the difference between the energy per atom {E) and the ground state energy per 
atom [Eq) in Hartree units vs. the number of integration steps. 
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Total energy minimization for a 64 atom liquid carbon sample, using steepest descent 
(SD), conjugate gradients (CG) and damped dynamics (D). The total energy corresponds 
to a parameterized tight-binding Hamiltonian (see text). We plot the logarithm of the 
difference between the energy per atom (E) and the ground state energy per atom (Eq) in 
Hartree units vs. the number of integration steps. The number of integration steps of the 
conjugate gradient calculation has been multiplied by two to take into account the increase 
in computational cost compared to the other methods. 

Total energy minimization for two randomized crystalline silicon samples using precon- 
ditioned steepest descent. The continuous line refers to a 32 atom cell with a long side of 
41 a.u. The dashed hne refers to a 64 atom cell with a long side of 82 a.u. We plot the 
logarithm of the difference between the energy per atom (E) and the ground state energy 
per atom (Eq) in Hartree units vs. the number of integration steps. In panel (a) we use a 
Yukawa potential to compute E and Eq, while in panel (b) the Yukawa potential is replaced 
by a Coulomb potential (see text). 
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